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Abstract 

We study the synchronizability and the synchronization dynamics of networks of nonlinear oscillators. 
We investigate how the synchronization of the network is influenced by some of its topological features 
such as variations of the power law exponent 7 and the degree correlation coefficient r. Using an ap- 
propriate construction algorithm based on clustering the network vertices in p classes according to their 
degrees, we construct networks with an assigned power law distribution but changing degree correlation 
properties. We find that the network synchronizability improves when the network becomes disassorta- 
tive, i.e. when nodes with low degree are more likely to be connected to nodes with higher degree. We 
consider the case of both weighed and unweighed networks. The analytical results reported in the paper 
are then confirmed by a set of numerical observations obtained on weighed and unweighed networks of 
nonlinear Rossler oscillators. Using a nonlinear optimization strategy we also show that negative degree 
correlation is an emerging property of networks when synchronizability is to be optimized. This suggests 
that negative degree correlation observed experimentally in a number of physical and biological networks 
might be motivated by their need to synchronize better. 

Keywords: Complex Network, Synchronization 

1 Introduction 

Networks of oscillators abound in physics, biology and engineering. Examples include communication net- 
works, sensor networks, neuronal connectivity networks, biological networks and food webs [Boccaletti, 
V.Latora, Y.Moreno, M.Chavez & D.U.Hwang 2006]. Under certain conditions such networks are known 
to synchronize on a common evolution, with all the oscillators exhibiting the same asymptotic trajectory. 
Synchronization was observed to play an important role in a wide variety of different problems (physical, 
ecological and physiological networks to name just a few); see for example [Y.Kuramoto 1984, A.T.Winfree 
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1980, G.B.Ermentrout & N.KopeU 1984, N.Kopell & G.B.Ermentrout 1986, X.F.Wang 2002, S.H.Strogatz 
2001, Mirollo & S.H.Strogatz 1990, S.H.Strogatz & Stewart 1993]. 

In this paper, we consider a network consisting of N identical oscillators coupled through the edges of 
the network itself [M.Barahona & L.M.Pecora 2002, T.Nishikawa, Motter, Y.Lai & F.C.Hoppensteadt 2003]. 
We suppose each oscillator is characterized by its own dynamics, x{t) = {xi{t),i = 1...N}, described by a 
nonlinear set of ODEs, x = f{x). The dynamics of each oscillator in the network is perturbed by the output 
function of its neighbors represented by another nonlinear term, say h = h{x). The equations of motion for 
each oscillator can then be given as follows: 



where a represents the overall strength of the coupling. Note that information about the network topology 
is entirely contained in the matrix C, whose entries Cij, j ^ i, are negative (zero) if node i is (not) connected 
to node j, with | Cij \ giving a measure of the strength of the interaction. 

In recent years, the analysis of large sets of data led to the identification of some important structural 
properties of many real- world networks; see for a review [Boccaletti et al. 2006]. Among these, the degree 
distribution P{k), with the degree k being the number of connections at a given node, was shown to be 
one of the most important features. In particular, scale-free networks, which are characterized by a power 
law degree distribution P{k) ~ k~^, have been observed to be widely spread in nature [Amaral, Scala, 
M.Berthlemy & H.E.Stanley 2000], [A.L.Barabasi & R.Albert 1999]. The main feature of scale free networks 
is an high heterogeneity in the degree distribution (higher than in purely random networks). Heuristically, 
this corresponds to the fact that real networks have often many low-degree nodes and only few nodes, termed 
as hubs, with many connections (thus leading to the high heterogeneity in the degree distribution). Such a 
characterization is particularly relevant since it has been shown to have important effects on the dynamics 
occurring over the network, such as the spreading of epidemics [Newman 20026], [R.Pastor-Satorras & 
A.Vespignani 2001] or the distribution of packets in traffic dynamics over the network [Goh, Kahng & 
D.Kim 2001], [M.Berthlemy 2004], [Arrowsmith, di Bernardo & F.Sorrentino 2005]. Also, it was shown that 
the scale-free structure of a network can have an important effect on its synchronizability (see Sec. 2 for 
further details). 

Recently, it has been suggested that real networks can often exhibit other important features. For in- 
stance, in the real world, nodes are not only abstract objects without any attribution, but each of them is 
characterized by some intrinsic features. Examples of interest could be the vertices age, spatial location, 
functional importance, level of activity and the number of connections each vertex has (that is the degree). 
An important topological property of physical and biological networks is that often their nodes show prefer- 
ential attachment to other nodes in the network according to their degree [Newman 2002 a], [Newman 2003]. 
Networks are said to exhibit assortative mixing (or positive correlation) if nodes of a given degree tend to 
be attached with higher likelihood to nodes with similar degree. (Similarly disassortative networks are those 
with nodes of higher degree more likely to be connected to nodes of lower degree.) 

The presence of correlation has been detected experimentally in many real-world networks. Interestingly, 
from the analysis of real networks, it was noticed that social networks are characterized by positive de- 
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gree correlation, while physical and biological networks show typically a disassortative structure [Newman 
2002a]. For example, in [A.Vasquez, R.Pastor-Satorras & A.Vespignani 2002], Internet was found to exhibit 
disassortative mixing at the Autonomous System level. 

In mathematical terms, degree correlation can be quantified by means of the observable r, introduced in 
[Newman 2002 a] and [Newman 2003]. In particular, the coefficient r was proposed in [Newman 2003] as a 
normalized measure of degree correlation, defined as the Pearson statistic: 

<Jq 

where qk is the probability that a randomly chosen edge is connected to a node having degree k, Uq is the 
standard deviation of the distribution qk and e,j represents the probability that two vertices at the endpoints 
of a generic edge have degree i and j respectively. 

The main aim of this paper is to investigate the relationship between the presence of degree correlation 
and the network synchronizability properties. Specifically, we shall seek to characterise how the presence of 
correlation on the network affects the synchronizability of the nonlinear oscillators at the nodes. We will 
show that correlation has indeed an effect on synchronizability in both weighed and unweighed networks. 
Our main finding is that disassortative networks synchronize better. Thus, as will be further discussed in the 
paper, the presence of negative degree correlation often detected in technological and particularly biological 
networks of nonlinear oscillators might be motivated by its benefits in terms of the network synchronizability. 

The rest of the paper is organised as follows. In Sec. 2 we briefly review the Master Stability Function 
approach and the effects on synchronizability of heterogeneity in the degree distribution, which is typical 
of scale free networks. In Sec. 3, we analyze how degree correlation can influence the synchronizability of 
unweighed networks and give an analytical explanation of the observed phenomena. Also, an optimization 
approach is used to motivate the conjecture that negative degree correlation can be an emerging property 
of networks when the aim is to optimize their synchronization. The case of weighed networks is discussed in 
Sec. 4, while the study of synchronization of both weighed and unweighed networks of Rossler oscillators is 
presented in Sec. 5. 

2 Effects of Topology on the Network Synchronizability 
2.1 The Master Stability Function Approach 

Recently, it has been proposed that the network topology, i.e. the way in which the oscillators are mutually 
coupled among themselves, has an important effect on the network synchronizability. This is defined in 
terms of the linear stability of the synchronous equilibrium {xi = X2 = ■■■ = Xn = Xg), whose existence is 
guaranteed by the fact that £ is a zero row-sum matrix. 

The stability of the synchronization manifold can be investigated by perturbing trajectories lying on it 
along directions which are orthogonal to the manifold itself. Namely, we consider sufficiently small pertur- 
bations Ci from the synchronous state so that Xi{t) = Xs{t) + ei{t), and: 
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de- ^ 

= Jf{xs)€i-a'^jCijJh{xs)€j, (3) 

where J denotes the Jacobian matrix. FoUowing the strategy presented in [L.M.Pecora & T.L.Carroll 1998], 
Eq. (3) can be transformed in blocks of the form: 

^ = [Jf{xs)-aXiJh{xs)]v, i = 1,2, ...,7V (4) 

where rj is defined in [L.M.Pecora & T.L.Carroll 1998] in terms of the perturbations e^, arc the real 
eigenvalues of the coupling matrix C, ordered in such a way that Ai < A2... < Ajv (in this manuscript we will 
deal only with real eigenvalues of the matrix C, although under general conditions, there may be complex 
conjugate eigenvalues). Note that Ai is structurally equal to 0, as the corresponding eigenvector is associated 
to the mode lying within the synchronization manifold. It is worth observing that using (4), the problem 
of studying the stability of a generic complex network in (1) has been replaced by that of considering the 
stability of N simpler independent systems. 

This approach is used in [L.M.Pecora & T.L.Carroll 1998] to derive the so-called Master Stability Function 
(MSF). Specifically, eq. (4) is considered as a parametric equation in a parameter a, given by: 

— = [Jf{xs) - aJh{xs)]v (5) 

and the values of a are sought corresponding to the maximum Lyapunov exponent of the system in (5) being 
negative. 

Interestingly, it can be shown that for a broad class of systems (associated to different dynamic functions 
/ and output functions h ), the Master Stability Function is negative in a bounded range of the parameter 

Of, say [oimini ^max] • 

Thus, in order to guarantee the stability of the synchronization manifold, all the crA,, for i = 2,3, ...,N, 
must lie in the range [cXmin,C(max]- Simply, this condition reduces to the following: 

aX2>amin 0-XN<amax, (6) 

or equivalently, to: 

Ajv 

^1 OLmin 

which guarantees the existence of at least one value of ct, for which the synchronization manifold is stable. 
Note that while OLmax and OLmin ai"© completely determined by assigning the functions / and h, Xn and 
A2 depend solely on the network topology. Thus, the synchronizability of a given network can be defined 
independently from the functional form of the dynamical systems at its nodes (i.e. of the functions / and 
h). 

In general, (6) defines a bounded range of values of a, say 1^ = [ "^J"- , ^f^], for which synchronization 
is attained. Therefore, an increase of A2 and a decrease of Ajv can lead to a larger interval I^- Thus, 
minimizing the eigenratio R = Xn/X2 yields a broadening of the range of values of a over which the network 
synchronizes. 
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Using the MSF approach, it is possible to investigate the effects on synchronization of the structural 
properties of the network topology (in terms of the eigenratio Hence, it is meaningful to characterise 

how the network topological features affect the Laplacian eigenratio. 

2.2 Synchronizability of Scale- free networks 

Sofar, the effects of the network topology on its synchronizability have been studied mainly with respect to 
the presence of scale free topologies or patterns in the network; see for example, [T.Nishikawa et al. 2003], 
[Motter, Zhou & J.Kurths 20056], [Fan & Wang 2005], [Lu, Chen & Cheng 2004], [Wu 2003], [Y.Moreno, 
Vazquez-Prada & F.Pacheco 2004], [Y.Moreno & F.Pacheco 2004], [Chavez, D.U.Huang, Amann, Hentschel 
& S.Boccaletti 2005, D.U.Huang, Chavez, Amann & S.Boccaletti 2005]. 

Scale free networks, which are common in nature, were found to show better synchronizability for increasing 
values of the power law exponent in [T.Nishikawa et al. 2003], [Motter et al. 20056]. Specifically, the rela- 
tionship was analyzed between the network structure and its synchronizability. An interesting phenomenon 
was observed which was termed as the "paradox of heterogeneity"; specifically, although heterogeneity in 
the degree distribution leads to a reduction in the average distance between nodes (the so called small world 
effect [D.J.Watts & S.H.Strogatz 1998]), it may suppress synchronization. 

To explain the observed phenomena, in [T.Nishikawa et al. 2003], the transition of the underlying network 
from scale free (power law distributed) to random (Poisson distributed) was shown to have a big impact 
on the eigenratio R of the Laplacian eigenvalues. Namely, a decrease of the heterogeneous nature of the 
network was discovered to yield, as a result, a reduction of R, thus increasing the synchronizability of the 
network itself. 

3 Effects of degree correlation on the Laplacian eigenratio 

In [di Bernardo, Garofalo & Sorrentino 2006], we proposed that a further decrease of R can be detected 
when negative degree correlation is introduced among the network nodes. Specifically, in [di Bernardo et 
al. 2006], using the configuration model [M.Molloy & B.Reed 1995], disassortative networks were found to 
synchronize better. In this section, we start by extending some of the results presented therein to the case of 
unweighed networks constructed by using the static model, which was firstly introduced in [Goh et al. 2001]. 

3.1 Network Construction methodology 

In this paper, in order to construct networks characterized by a given degree distribution, we will use the 
following algorithm: 

1. We associate to each vertex i = a weight or fitness (jji oc {i + 6)~", with a = [0, 1), 6 being a 
fixed parameter. 

2. Then, links are added among the network nodes with probability proportional to the vertices weights, 
until £ links have been introduced. 
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Figure 1: N / J\f as varying the degree distribution exponent 7 in networks characterized by different values of 
the degree correlation coefficient r. The legend is as follows: r = —0.3 (o), r ~ O(-), r = +0.3(A). M — 10^ 
nodes, f = 4 ■ 10^, 9 = 10. 

Note that by using this methodology, the expected degree at each node i, is ki ~ (piE^ and it can be shown 
that the degree distribution is expected to follow a power law, P{k) ^ k with 7 = 

Observe that the nodes fitness (pi , is the only parameter of the model; we wish to emphasize that assigning 
the fitness corresponds to assigning the expected degree at the vertices, thus the degree at the nodes can be 
considered as an equivalent parameter of the model. 

By using such a model, the giant component of the network may not include all the network nodes, i.e 
the number of connected nodes, say iV, may be < M. In Fig. 1, N/Af has been plotted as varying both 7 
and r. In what follows, we will consider only the largest connected part of each generated network and refer 
to N as the number of nodes belonging to the giant component. 

Using this methodology to construct a network, we can handle the transition from very heterogenous scale 
free networks (when a — 1 and 7 = 2), to highly homogenous ones (as a approaches and 7 approaches 
cxi); note that in the limit of 7 ^ 00, in which to each vertex is associated the same weight, we recover the 
classical random network model, which was introduced in [Erdos & A.Renyi I960]. Moreover, this does not 
imply a dependence of the total number of edges on the degree distribution exponent 7 (as is the case of the 
widely used configuration model [M.MoUoy & B.Reed 1995]). 

By using this model we are able to reproduce networks characterized by: (i) a given number of nodes TV, 
(ii) a desired number of edges £, (iii) a degree distribution with a desired value of the exponent 7 (note 
that the same model has been already employed to study the synchronization properties of networks in [Lee 
2005] and in [Oh, Rho, Hong & B.Kahng 2005].) Then it is possible to devise a strategy similar to the one 
presented in [Newman 2003], to generate networks with a given level of assortativity, i.e. a desired value of 
the degree correlation coefficient r. 

Notice that the model introduced in [Goh et al. 2001] is known to exhibit spontaneous negative degree 
correlation [Lee, Goh, Kahng & Kim 2006] when 2 < 7 < 3 (this phenomenon being referred to in [Lee et 
al. 2006] as a fermionic constraint, due to the impossibility of connecting pairs of nodes with more than 
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one link). The same phenomenon has been observed also in networks constructed by using the configuration 
model, as described in [J.Park & M.E.J. Newman 2003, D.S.Callaway, J.E.Hopcroft, J.M.Kleinberg, Newman 
& S.H.Strogatz 2001, Catanzaro, Boguna, & Pastor-Satorras 2005]. However, in what follows, we will 
disregard this particular effect; namely, we will compare networks characterized by different degree correlation 
properties according to the values of the index r (including the case r = 0), which have indeed been measured 
over the networks. 

3.2 Unweighed networks 

In the rest of this section we will investigate the case of unweighed network topologies, with all the links 
being associated to a constant unitary weight. Specifically in such a case, the matrix £ can be written as 
follows: 

jC = D-A, (8) 

where D — {-Dy } is a diagonal matrix such that Da = ki, i = 1, 2, N, and A is the associated adjacency 
matrix. The case of weighed topologies will be considered in Sec. 4. 

The overall results arc summarized in Fig. 2(a) where the effects of varying the degree correlation on the 
Laplacian cigcnratio R arc shown for different values of the degree distribution exponent 7. As discussed 
in [T.Nishikawa ct al. 2003], in the case of uncorrclatcd networks (r = 0). synchronizability improves for 
increasing values of 7. Moreover, for all values of 7, we observe a reduction of R for decreasing values of r. 
This means that disassortative mixing enhances the network synchronizability. Interestingly, as depicted in 
Fig. 2(b) and Fig. 2(c), we observe that, under variations of the correlation parameter, the changes in R 
seem to be mainly due to variations of A2. This is not surprising if we consider that A^v is known to scale 
with kmax [B.Mohar 1989o], the maximum degree at the vertices and kmax cannot vary with the degree 
correlation coefficient r. 

It is worth noting here that the results provided in this section are obtained by using a different network 
generation model from that used in [di Bernardo et al. 2006] . Moreover, the results shown in this paper (and 
in particular those in Fig. 2) are found to confirm those previously presented therein. 

3.3 An analytical explanation 

Here we briefly expound in this context the analysis provided in [di Bernardo et al. 2006], in order to provide 
an explanation of the behavior of the eigenratio i? as a function of the degree correlation coefficient r depicted 
in Fig. 2. 

Suppose the network vertices are divided in p classes according to their degree, i.e. clustering in each 
class C, all the vertices having degree fc^, (i — 1, 2, ...p). Let us term as Ui the number of vertices in class Cj, 
then the probability of finding a vertex belonging to class Ci at the end of a randomly chosen edge within 
the network is given by Qi = riiki/ ^^n^k^. In so doing, following [Newman 2002a], the presence of degree 



7 



correlation can be estimated by using the coefRcient r defined as: 

_ k^(E-qq^)k 



(9) 



where Ug is the standard deviation of the distribution g;, k — (fci, fc2, /cp)'^, q = ((71, 52, 9^)"^ and 
E — {cy} G E^'^P, with Cij being the probability that a randomly chosen edge in the network connects 
nodes having degree i and j . 

From (9), it is possible to obtain the distribution of edges among the network vertices as a function of r 
as follows: 

E = qq"^ + rff^M, (10) 

where M is a symmetric matrix having all row sums equal to zero and appropriately normalized such that 
k'^Mk = 1. Specifically, we can express M as follows: 



M 



mm 



(k^m)2 ' 

where m ~ (mi, m2, mp)^ is a vector such that X^i^'^i — i^'^^ instance we can choose m such that 
TTii < mj_|_i for j = 1, 2...,p— 1 in order to have a convenient form of the matrix M with positive values near 
the main diagonal and negative values far away from it). 
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Figure 2: Synchronizability of degree correlated scale free networks of size 10"^ nodes. Behavior of the eigenratio 
A]v/A2 (a), of the second lowest eigenvalue A2 (b) and of the highest one Ajv (c), as functions of the correlation 
coefficient r defined in [11], for 7 varying from 2.4 (top line) to 4.2 (bottom line) in steps of 0.3. (d) The eigenratio 
as function of 7 as varying the correlation coefficient r from —0.3 (bottom line) to 0.3 (top line) in steps of 0.1. The 
lines are guided by the eye. 
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Figure 3: The isoperimetric problem consists in finding the subset S such as to minimize the ratio hoiS) in 
Eq. (11). In this schematic picture, a given subset S is selected as the set of all network vertices contained 
inside the circle in the figure. Its cardinality is given by the number of nodes in the circle, \S\ = 5 and its 
boundary is given by the number of edges (represented in bold) which connect the subset S to the rest of 
the network, V{S) = 9. 

As Ajv is fomid to bo almost iuclcpcudent from the correlation coefficient r (see Fig. 1(c)), we will focus now 
on estimating the effects of correlation on the eigenvalue A2, the parameter known as algebraic connectivity 
of graphs [B.Mohar 1989&]. Specifically, to estimate an upper bound for A2, we use the following Cheeger 
inequality from graph theory [J. Cheeger 1970]: 

A2 < ho, (11) 

with ho, the Cheeger constant of a graph, defined as follows [J. Cheeger 1970]: 

ho =mmhG{S). (12) 

For a given subset of vertices, say S, hoiS) is the quantity given by: 

V{S)N 

'^^^^^ - |51(Ar-]5|)' (^^^ 
where ©(5) is the number of edges in the boundary of S and 15] < ^, is the number of vertices in S. 

The problem of finding ho, known as the isoperimetric problem in graph theory, is illustrated by means 
of a representative example in Fig. 3. Note that, given a graph of order N, there are C{N, i) possible 

ways of choosing the subset S, where C{a, b) is the number of combinations of a elements in b places, and 
thus the problem of choosing S such as to minimize haiS) is NP-hard. 

We will show that (11) can be successfully used to compute an upper bound on A2. To overcome the 
limitations due to the computation of the subset S that minimizes haiS), we will follow a stochastic approach 
in order to estimate hciS), starting from the available information we have on the network. We will assume 
that the noticeable features of the network are only the degree distribution and the correlation specified; all 
other aspects being completely random. 

Note that, by considering as equivalent all the vertices characterized by the same degree, the original 
problem in Fig. 3, can be converted in the other (simpler) one shown in Fig. 4. In the picture, the subset S 
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Figure 4: The problem described in Fig. 2 has been converted in the one represented in this picture. Here, 
the subset S is identified by the degrees of the vertices inside it, which are two nodes of degree 3, two nodes 
of degree 4 and one node of degree 7. Note that, in the original formulation, the problem was of choosing 
the vertices to be included in S in order to minimize hoiS). Here we have restated it by considering only 
the degrees of the vertices to be included in S, in terms of the vector y, such as to minimize hciy)- 

includes two vertices of degree 3, two vertices of degree 4 and one vertex of degree 7. In terms of the nodes 
degrees the original problem, consisting of finding the optimal set of vertices, can then be re-interpreted as 
that of finding the optimal mix of degrees of nodes in S. 

Then, we can give a full characterization of a randomly chosen subset S in terms of the number of nodes 
in it, say Xi, belonging to each class Ci {i = 1,2..., p). Let us term as yi = ^ the fraction of nodes in 
S drawn from each class Ci (i = 1,2..., p). Then problem (12) becomes that of finding the combination 
y = {yi,y2, ■.■,yp)^ such as to minimize hoiy), i.e. miuy /iG(y, r). Then, since under this formulation 
nodes with the same degree arc all clustered in the same class, the complexity of the problem is reduced 
to approximately "^fJi C{N,i)/Il'^!!^^''nj. It is also worth noting here that the subset S is not supposed to 
satisfy any particular condition, not even of being connected. 

Now, we observe that the number of edges in the boundary, say ^{S), is given by the total number of 
edges starting from the vertices in S, less the ones, say 1{S), that are contained in S, i.e. having both 
endpoints in S. Thus we can estimate I and V as follows: 

I{S)=I{y,r)=£y^Ey, 
V{S) = V{y, r) = x^k - 2J(5) = 2£(y^q - y^£;y), 

where x = {xi,X2, Xp)'^ and £ is the total number of edges in the network. 
Therefore hciS) becomes: 

, , , 2£(y'^(i-y'^Ey)N 

under the constraint that n-^y < N/2, where n is the vector (m, n2, rip)^. A numerical optimization 
algorithm can then be used to find the subset S that minimizes hG{y,r) in terms of yi,y2, ■■■,yp (and 
subsequently r) and, in turns, an upper bound for A2. 
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Also, from (10) and (14), we get: 

dhciS) dV(S) „^ 2. T ^2 „ 

oc = -2£a2(y^m) < 0. (15) 

Since, for any vector y, (15) is satisfied, then we have that < 0. Therefore, we can predict analytically 
that he and hence A2 will be decreasing as the degree correlation is increased and, as a consequence, the 
eigenratio ^ will increase for higher values of the correlation coefficient. 

Another interesting inequality in spectral geometry is due to Mohar [B.Mohar 1989a]: 



A2 > kmax ~ Y kmax^ ~ ^'q , (16) 

where h'^ = ming = 2£(y^q — y'^ Ey) / [n^ y) . Using (16), we can also get a lower bound on A2. Then 
following an approach similar to the one used to compute the upper bound, it is easy to show that the lower 
bound in (16) has to decrease with r (note that when making the correlation change, the degree distribution 
is fixed and thus, kmax cannot vary with r). Since both the upper and the lower bounds have to decrease 
with r, A2 is also expected to have the same trend. 

3.4 Disassortativity as an emerging property 

The main result of the derivation presented above is the finding that disassortative networks synchronize 
better. We wish now to assess whether negative degree correlation can be thought of as an emerging property 
of networks with an assigned degree distribution in order to improve their synchronization. In particular, 
we noticed that the variation of the degree correlation affects mainly the Laplacian eigenvalue A2. Thus, we 
shall seek to find if varying the correlation, while keeping the degree distribution fixed, is indeed a good way 
of optimizing the network synchronizability properties. 

This might be a solution to a classical problem in graph theory optimization which is the construction of 
expander graphs, i.e. highly efficient communication networks, characterized by high values of A2 [N.Alon 
1986] (and therefore low values of the eigenratio R). 

We will use a simulated annealing meta-heuristic technique, to solve the problem of maximizing A2 while 
keeping unchanged the network degree distribution. To this aim, given a network with a certain degree 
distribution, we will perform the following iterative procedure. At each step, the endpoints of a randomly 
selected pair of edges are exchanged if exp 

(^MM) > 2, where z is an uniformly distributed random variable 
between and 1, A(A2) is the variation achieved in the objective function A2 before and after the execution 
of the move and T is a control parameter, which similarly to the original formulation of the simulated 
annealing procedure is known as the system 'temperature'. As the algorithm runs, the temperature T is 
decreased according to an exponential cooling scheme (see [S.Kirkpatrick, C.D.J.Gerlatt & M.P.Vecchi 1982] 
for further details). 

As shown in Fig. 5, while trying to maximize A2, we observe the spontaneous emergence in the network 
of interest of negative degree correlation, namely an increase of its disassortativity. 

Thus, following an entirely different approach, we come to the conclusion that if the degree distribution 
of a given network is fixed, then to improve its synchronizability one has to introduce negative degree 
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Figure 5: Effects of the maximization of A2 on the correlation parameter r defined in [11] for different scale- free 
networks with 7 = [2, 3, 4, 5]. 



correlation among its node. This suggests that in evolutionary biological networks of nonlinear oscillators, 
disassortativity might be an emerging property necessary to optimize the synchronization process. 



4 Weighed Networks 

In this section we shall extend the results presented above to the case of weighed topologies. Interestingly in 
[Motter et al. 20056, Chavez et al. 2005], it has been shown that an appropriate choice of weights over the 
network links, may improve considerably the network synchronizability. In [Chavez et al. 2005], it has been 
shown that such improvement is highest if the weights are set to be proportional to a particular topological 
property at the network links, termed as load or betweenness centrality. 

We have seen sofar that scale-free networks can synchronize better as the degree distribution exponent 
7 is increased and the degree correlation coefficient r is decreased. In [Motter, Zhou & J.Kurths 2005a], it 
was shown that synchronizability of such networks can be further enhanced by an appropriate choice of the 
network coupling weights, that takes into account exclusively topological information. Specifically, it was 
found that in order to enhance the synchronizability of heterogeneous networks, the coupling matrix should 
have the form: 

C = D-f^iD - A), (17) 

where the strength of the coupling, for each edge, is scaled by a power (with exponent f3) of the degree of 
the starting node. Note that, by tuning the parameter /3, it is possible to vary the strength of the coupling 
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p p p 

Figure 6: Plots of \n/^2 as varying /3 for scale free networks with different values of the power-law exponent 
7 — [2.5, 3, 3.5]. On the left plot we show negatively correlated networks (r = —0.3), on the right positively correlated 
ones (r — +0.3), in the center, networks characterized by absence of degree correlation (r = 0). The thick lines 
represent the behavior of the eigenratio for random networks (7 = 00). A'' = 10^ nodes, £ — 4 ■ 10^ . 

from high to low degree vertices and viceversa. Also, from (17), C can be rewritten as D~^/^{D — A)D~^/'^ 
and thus, even if £ is not symmetric, its spectrum is always real, whatever the value of /3. 

In [Hotter et al. 2005a], in the case of heterogeneous networks, an optimal coupling was found to be /3 = 1 
in (17). It is worth noting that this leads to asymmetric coupling. Take for example the link between an 
high degree and a low degree vertex; a coupling of the form (17), with /3 = 1, results in the low degree node 
having a much higher influence on the higher degree one than viceversa. 

Generally speaking, the case oi (3 — 1, with the weights of the incoming links at each vertex summing to 
-1 (so that Cii ~ 1 Vi), is of particular interest in that all the networks satisfying this constraint, become 
directly comparable among themselves, in terms of their spectral properties [Chavez et al. 2005] (see also 
[Sorrentino, di Bernardo, Cuellar & Boccaletti 2006].) In fact, according to the Gerschgorin's circle theorem, 
all the eigenvalues (Aj — X^+jXj, j — 1, 2, N) of these network are known to belong to a bounded region 
of the complex plane, and specifically the circle of radius 1, centered at 1 on the real axis [Chavez et al. 2005]. 
Thus all the A^, are constrained to belong to the interval [0,2] (i.e. = X[ < X2 < ... < A^ < 2) and the A* 
to the interval [-1,1], j = 1, 2, iV (in the case of interest here of eq. (17), A* = Vj). Moreover the case 
of /3 = 1 has been shown to represent an optimum in terms of the network synchronizability, independently 
from the form of its degree distribution. 

Moreover in [Motter et al. 2005a], it was claimed that in the optimally coupled network (/? = 1), synchro- 
nizability is solely determined by the average degree, while heterogeneity in the network connectivity does 
not affect synchronization. However, as shown in Sec. 4.1, we observe that in the case of assortative mixing, 
this is not true anymore; specifically in the particular case of r = 0.3 shown in the right panel of Fig. 6, 
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Figure 7: Synchronizability of degree correlated networks in the optimal regime {(3 — 1). Left figure: behavior of the 
eigenratio Ajv/A2, as a function of 7. N = 10^, 5 = 4- 10"^. Legend is as follows: r = — 0.3(o), r = — O.l(x), r = 0(*), 
r = O.l(V), r — 0.3(0)- Right side: behavior of the eigenratio Ajv/A2 as a function of the correlation coefficient r, 
for 7 varying from 2 to 5. A'' = 10^, £ = 4 - 10^. Legend is as follows; 7 = 2(-) , 7 = 3(+), 7 = 4(*), 7 = 5(x). 

homogenous networks are found to behave better than heterogenous ones, even in the optimal case oi P — 1. 
4.1 Effects of degree correlation 

Here, we will show that, assuming a coupling of the form (17), the presence of degree correlation does not 
alter the optimal value for (3 and better synchronizability properties are still detected when (3 — 1 as was 
observed in the uncorrelated case. Moreover, the presence of degree correlation, unlike heterogeneity, does 
continue to have an effect on the network synchronizability, even in the optimal regime. (Namely, as in the 
case of unweighed coupling, disassortative mixing always enhances network synchronization). 

In Fig. 6 we have reported the behavior of the eigenratio as varying /3 in networks with different degree 
correlation properties (characterized by r = [— 0.3, 0, +0.3]). Our numerics confirm that the eigenratio 
R = Ajv/A2 is characterized by a peaked minimum at /3 ~ 1, independently from the value of r (i.e. also 
when r ^ 0). Note that in the case of uncorrelated networks, when /3 = 1, the eigenratio is quite insensible 
to the specific form of the degree distribution, as claimed i.e. in [Motter et al. 2005a]. Nonetheless, as r 
increases, the eigenratio appears to be more sensible to variations in the degree distributions, even in the 
optimal case of /? = 1. Specifically, as shown in Fig. 7, in such a case the minimum is more pronounced for 
higher values of 7 (i.e. the paradox of heterogeneity is still present even at /? = 1). 

Differently from the case of unweighed networks, we notice that in Fig. 8 Xn is now sensible to variations 
in r, when compared to A2; thus the behavior of the eigenratio R as r varies, cannot be explained only 
in terms of the second smallest eigenvalue but is a combination of variations of both the eigenvalues (see 
Fig. 8(a) and 8(c)). The behavior of A2 and Xn under variations of 7 for different values of r is shown in 
Fig. 8(b)-(d). Namely, when increasing 7, we observe a rise in A2 and at the same time, a decrease in Ajv, 
leading to an overall decrease of the eigenratio R and therefore better synchronizability. Similarly, when r 
is varied from the disassortative (r < 0) to the assortative (r > 0) regime, we observe both a decrease in A2 
and an increase in Aat, leading to a rise in the eigenratio R. 
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Figure 8: Synchronizability of degree correlated networks in the optimal regime (/3 = 1). N = 10 nodes, £ — 4 - 10 . 
Behavior of the highest eigenvalue Ajv (a) and of the second lowest one A2 (c), as functions of the correlation coefficient 
r, for 7 = [2,3,4,5]. Behavior of the highest eigenvalue Xm (b) and of the second lowest one A2 (d), as functions 
of 7, for r = [-0.3,0, +0.3]. Legend in (a) and (c) is as follows: left side, r = -0.3(o), r = -0.2(+), r = -O.l(x), 
r = 0(*), r = O.l(V), r = 0.2(-), r = 0.3(0). Legend in (b) and (d) is as follows: 7 = 2(-) , 7 = 3(+), 7 = 4(*), 
7 = 5(x). 

Observe that in the case of the normaUzed Laplacian (17) considered here, the following relationship from 
graph theory is valid [Chung 1997]: 



< A2 < 2h"G, 



(18) 



where h'^ = mins 11^(3) and hQ{S) is defined as: 



V{S) 



min {vol{S), vol(S)) 
with vol{S) — J2ieS ^ote that in this case hQ{S) can be estimated as: 



(19) 



h"G{y.r) = — 



min(y^q, (1 - y'^)(\) ' ^"^^^ 
Then it is sufficient to observe that, for each S, the denominator in (20) is a quantity independent from 
r, and it becomes straightforward to obtain the existence of an upper and a lower bound for A2, decreasing 
with r. Unfortunately notwithstanding this, estimating the effects of r on Aat resulted to be a non-trivial 
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task and therefore it remains an open problem. An important feature seems to be the symmetry between 
the behavior of A2 and Ajv reported in Fig. 8(a)-(c) and 8(b)-(d). 

From the numerics shown in Fig. 8, we observe a clear advantage of introducing negative degree correlation 
within the network (in terms of both Xj^ and A2), when /3 = 1. Thus we expect to observe an enlargement of 
the synchronization interval I„, for both higher and lower values of a as decreasing r. Moreover, since this 
happens for whatever a choice of /?, we find that disassortative mixing always enhances synchronization. 

The effects of correlation on the eigenratio R in the case of /3 = 1 (i.e. optimal synchronizability) is fully 
illustrated in Figs. 7 and 8. Specifically from Figs. 7 and 8, we observe that negatively correlated networks 
continue to show higher synchronizability, for every value of the exponent 7. This indicates that the effects 
of correlation on R are not dampened by the presence of an appropriate coupling of the form (17). Fig. 
7 also shows a non-negligible dependence of the eigenratio on the exponent of the power-law distribution 
7, indicating that heterogeneity continues to affect the synchronizability of weighed correlated networks, 
especially in the cases where r > (as also shown in the right panel of Fig. 6). 

5 Effects of degree correlation on the synchronization dynamics 

Now, wc shall seek to investigate the effects of variable degree correlation on the dynamics of networks of 
oscillators. Specifically in what follows we will evaluate the behavior of dynamical networks of (i) identical 
and (ii) non identical oscillators. 

As an example of identical oscillators, we consider a network of coupled chaotic Rossler oscillators, con- 
nected as in (1). In [di Bernardo et al. 2006], we have already reported numerical results showing the 
synchronization dynamics of such networks in the case where h{x) — Ix, with / being the 3-dimensional 
identity matrix. Therein we observed that negative degree correlation does indeed have beneficial effects on 
the network synchronization. Specifically, our preliminary results confirmed that negative degree correlation 
is able to enhance the network synchronization dynamics. 

Here we suppose the dynamics at each node i to be described by the following vector field, Vi = {xi,yi, Zi): 



where a is the tunable strength of the coupling. Note that, according to this scheme, we have chosen the 
Rossler to be coupled through the variables x and z. The reason is that under this assumptions, the MSF 
results to be negative in a bounded range of values of the parameter a, yielding that the coupling a should 
be neither too low, nor too high. 

We have performed simulations over scale-free complex networks, with variable power-law exponents, 
characterized by several values of the degree correlation coefficient (typically in the range r G [—0.3, -1-0.3]). 
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Figure 9: Synchronization dynamics over unweighed networks (/3 = 0). The plots show the order parameter 
at regime m for networks characterized by different 7 (in each subplot) and r (differently colored lines) for 
increasing values of the coupling gain cr; N = 10^, f = 4 • 10^. Note the different scale on the y-axis in the 
lower figure. 

In order to have a measure of the overall network synchronization, we introduce the order parameter 
mit) e [0, 1], defined as follows (see also [Yook & Meyer-Ortmanns 2005]): 



where Q{x) is the Heavyside function, i.e. Q{x) = 1 if x > and 0(x) = otherwise. The parameter 5 is 
a small number to account for the finite numerical accuracy, (here we have chosen S — 10~^), so that two 
states in phase space lying inside a sphere of radius 8 are considered as mutually being synchronized. The 
parameter m{£) gives the fraction of pairs of elements (i, j) which are synchronized at time t (i.e., dij < S). 
This fraction is equal to unity if all possible pairs are synchronized and zero if no pair is synchronized, 
with intermediate values < r < 1 indicating partial synchronization. In what follows we will look at the 
asymptotic average value of m{t), say m, as function of the coupling strength cr. 

Specifically, in Figs. 9 and 10, a double phase transition is observed: the first one from a non-synchronous 
to a synchronous regime as the coupling strength is increased; the second one from the synchronous to a 
non-synchronous phase as the coupling strength is further increased and exceeds a critical value. 

At first, we suppose the network topology to be unweighed, i.e. the coupling is identical over all the 
edges and equal to 1 . In order to explore the effects of the network structure (degree distribution and degree 
correlation) on the dynamical synchronization process, we have performed several numerical simulations, the 




(21) 
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Figure 10: Synchronization dynamics over weighed networks {f3 ~ 1). The plots show the order parameter 
at regime m for networks characterized by different 7 (in each subplot) and r (differently colored lines) 
for increasing values of the coupling gain a; N ^ 10"^, £ — 4 ■ 10'^. The lower plots are enlargements of 
the upper ones. In the right-low panel, it is interesting to observe how the range of values of a for which 
synchronization is achieved (at m = 1) increases as r is reduced from 0.3 to -0.3. 

main results being illustrated in Fig. 9. The asymptotic value of the order parameter m, has been computed 
for networks with different topological features. We have compared scale free networks characterized by 
different power-law exponents and we found that, as can be observed by comparing the plots in Fig. 9, 
homogenous networks (characterized by higher values of 7) synchronize better than heterogenous ones (note 
the different scale on the y-axis in the lower plot, where heterogeneous scale-free networks, with 7 = 2 
cannot be synchronized in practice for any value of a). This indicates that an high heterogeneity in the 
degree distribution can result in being very harmful to the overall synchronization [T.Nishikawa et al. 2003]. 
Moreover negative degree correlation is seen to enhance the network synchronizability, for all the values of 
7 considered (excluding those for which synchronization cannot be achieved, e.g. 7 = 2 in Fig. 9, where 
TO ^ 1). This is consistent with the predictions based on the eigenvalue analysis provided in Sec. 3 and 
confirms the beneficial role played by negative degree correlation with respect to networks synchronization 
even in terms of synchronization dynamics. 

The case of weighed topologies, (i.e. when /3 = 1) is shown in Fig. 10. Again, negative degree correlation 
is observed to enhance the network synchronization, at high values of the order parameter when m « 1. 
Moreover, as shown in Fig. 10, the network synchronization is reduced as the heterogeneity in the degree 
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distribution is increased (i.e for lower values of the exponent 7). This indicates that, even in the optimal 
regime where (3=1, networks characterized by variable degree distributions, behave differently in terms 
of their synchronization dynamics (practically we observe an analogous phenomenon as in the case of the 
unweighed topologies). Also, it is interesting to observe in Fig. 10 that the range I„ of values of a for which 
synchronization is completely achieved (at m = 1), increases by reducing r from 0.3 to -0.3. 

It is worth noting here that, by comparing Figs. 9 and 10, weighed networks, in which the strength of 
the coupling has been rescaled by the degree at each node, are more synchronizable than unweighed ones 
(as also predicted by the behavior of the eigenratio in Fig. 6). We wish to emphasize that all the numerical 
results provided in this section, are in good agreement with the eigenvalue analysis presented in Sec. 3. 

6 Conclusions 

The structure of many real world networks is characterized by non-trivial degree-degree correlation in the 
network connections. This leads alternatively to disassortative mixing in biological and technological net- 
works and assortative mixing in social networks. In this paper wc have studied the effect of degree correlation 
on the network synchronizability and synchronization dynamics. Specifically, we studied the effects of corre- 
lation on the Laplacian eigenratio, a parameter proposed in [M.Barahona & L.M.Pecora 2002] as a measure 
of the synchronizability of a network of coupled nonlinear oscillators. 

We have found that disassortative mixing, which is typical of biological and technological networks, plays 
a positive role in enhancing network synchronizability. The numerical observations were confirmed by the 
analytical estimates found for the Laplacian eigenratio, which were shown to be well suited to describe the 
observed phenomena. 

Following [Motter et al. 2005 a], we then analyzed the effects of weighted and directed coupling of the 
form (17) and found that the presence of correlation continues to affect synchronization with disassortative 
weighted networks synchronizing better than assortative ones. We noticed that, even in the presence of 
degree correlation, synchronizability seems to be optimal when the strength of the coupling along each edge 
is made inversely proportional to the degree of the starting node. Though a common belief is that when 
such coupling is considered, heterogeneity becomes unable to siippress synchronization, our numerics show 
that an exception to this paradigm is represented by the case of assortatively mixed networks, where the 
effect of the degree distribution is found to be strongly enhanced. 

Using a nonlinear optimization approach we found that negative degree correlation is naturally attained by 
the network when the aim is to minimize A2 and hence enhance^ synchronizability. Thus, we conjectured that 
disassortative mixing has played the role of a self organizing principle in leading the formation of many real 
world networks as the Internet, the World Wide Web, proteins interactions, neural and metabolic networks. 

Finally, we investigated the synchronization dynamics of both weighed and unweighed networks of identical 
Rossler oscillators, confirming the theoretical results obtained in the paper. A more general case of weighed 
degree correlated networks (with complex spectrum of the Laplacian) will be discussed elsewhere [Sorrentino 
et al. 2006]. 
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